function [ v, zeta , Lh , L , vL , C , ce , ...
           pop , f , w , wh , u , H , h , s ] = ...
    SolveEntry( Grids , Params , GE )

% This function solves for entry using bisection
% unemployment rate is decreasing in entry cost

ScaleFactor = Grids.ScaleFactorEntry ;
run MiscFunctions.m

%% Initialize bounds

ce0 = Params.ce ;
ce1 = Params.ce ;

%% Lower bound

disp(' ')
disp('************')
disp('Looking for lower bound of entry cost bisection...')
disp('************')
disp(' ')

Params.ce = ce0 ;

[ v , zeta , Lh , L , vL , C ] = ...
    SolveC( Grids , Params , GE ) ;

% update C for speed
GE.C = C ;

MeanU = uMean(v,zeta,L) ;

ceratio0 = MeanU / Params.TargetUrate ;
enteredloop0 = 0 ;

% unemployment rate is an increasing function of the entry cost

if abs( ceratio0 - 1 ) > Grids.tolu
    
    while ceratio0 > 1

        disp(' ')
        disp('************')
        disp(['ceratio0 = ' num2str(ceratio0)])
        disp(['ce0      = ' num2str(ce0)])
        disp('************')
        disp(' ')

        ce0 = ce0 / ScaleFactor ;
        Params.ce = ce0 ;

        [ v , zeta , Lh , L , vL , C ] = ...
            SolveC( Grids , Params , GE ) ;

        % update C for speed
        GE.C = C ;

        MeanU = uMean(v,zeta,L) ;
        ceratio0 = MeanU / Params.TargetUrate ;

        enteredloop0 = 1 ;

    end

    %% Upper bound

    disp(' ')
    disp('************')
    disp('Found lower bound of entry cost bisection')
    disp('Looking for upper bound of entry cost bisection...')
    disp('************')
    disp(' ')

    if enteredloop0 == 1
        ce1 = ScaleFactor * ce0 ;
        Params.ce = ce1 ;
    else
        Params.ce = ce1 ;
    end

    [ v , zeta , Lh , L , vL , C ] = ...
           SolveC( Grids , Params , GE ) ;

    % update C for speed
    GE.C = C ;

    MeanU = uMean(v,zeta,L) ;
    ceratio1 = MeanU / Params.TargetUrate ;


    enteredloop1 = 0 ;      
    while ceratio1 < 1

        disp(' ')
        disp('************')
        disp(['ceratio1 = ' num2str(ceratio1)])
        disp(['ce1      = ' num2str(ce1)])
        disp('************')
        disp(' ')  

        ce1 = ScaleFactor * ce1 ;
        Params.ce = ce1 ;

        [ v , zeta , Lh , L , vL , C ] = ...
           SolveC( Grids , Params , GE ) ;

        % update C for speed
        GE.C = C ;

        MeanU = uMean(v,zeta,L) ;
        ceratio1 = MeanU / Params.TargetUrate ;

        enteredloop1 = 1 ;

    end

    if enteredloop1 == 1
        ce0 = ce1 / ScaleFactor  ;
    end

    %% Bisect

    disp(' ')
    disp('************')
    disp('Found both initial values for bisection of entry cost')
    disp(['ce0      = ' num2str(ce0)])
    disp(['ce1      = ' num2str(ce1)])
    disp(['ceratio0 = ' num2str(ceratio0)])
    disp(['ceratio1 = ' num2str(ceratio1)])
    disp('************')
    disp(' ')

    %pause

    error = 1 ;
    it = 1 ; 
    while error > Grids.tolu && it < Grids.maxit

        cemid = ( ce0 + ce1 ) / 2 ;
        Params.ce = cemid ;

        [ v , zeta , Lh , L , vL , C ] = ...
               SolveC( Grids , Params , GE ) ;

        MeanU = uMean(v,zeta,L) ;
        ceratio = MeanU / Params.TargetUrate ;

        error = abs( ceratio - 1 ) ;

        disp(' ')
        disp('************')
        disp('Bisecting on entry cost')
        disp(['ce      = ' num2str(cemid)])
        disp(['ceratio = ' num2str(ceratio) ])
        disp('************')
        disp(' ')


        if ceratio > 1
            ce1 = cemid ;
            ceratio1 = ceratio ;
        elseif ceratio < 1
            ce0 = cemid ;
            ceratio0 = ceratio ;
        end

        it = it + 1 ;
        %disp([ 'it = ' num2str(it) ] )

    end

    frac = ( ceratio - ceratio0 ) / ( ceratio1 - ceratio0 ) ;
    ce = ( 1 - frac ) * ce0 + frac * ce1 ;
       
    Params.ce = ce ;

    [ v , zeta , Lh , L , vL , C ] = ...
           SolveC( Grids , Params , GE ) ;
        
else
    
    ce = Params.ce ;

end
 

%% CONSTRUCT OUTPUT

% Recover actual population density
L = Lh / sum( Lh .* Grids.mx ) ;

pop = L .* Grids.mx ; pop = pop / sum(pop) ;

x0 = Params.muLearn / ( Params.muLearn + Params.DeltaLearn ) ;                

run MiscFunctions.m

f  = Phi(v,zeta) ;          
w  = NomWagesNoHC(v,zeta) ;
wh = NomWagesHC(v,zeta) ;
u = uRate(v,zeta) ;
H = MeanHC(x0,Lh,v,zeta) ;
h = Params.DLearn ./ ( Params.DLearn + Params.phi * u ) ;
s = Params.delta * zeta ;




end

